#include <cstdio>
#include <cmath>
#include <vector>
#include <queue>
#include <algorithm>
#include <bitmap_image.hpp>

unsigned int Height, Width, Size;

const unsigned int CNTNUM = 32;
unsigned int CNTWIDTH = 256 / CNTNUM;

const unsigned int LINEMUL = 65536;

struct point {
	float Sum, Sumx, Sumy, Sumxy, Sumx2;
	float Avgx, Avgy, Avgxy, Avgx2, Ver, Par, Diffy, width, Sig;
};

int cnt[CNTNUM];
std::vector <unsigned int> black;
std::vector <bool> light;
float Sum, Sumx, Sumy, Sumxy, Sumx2, Avgx, Avgy, Avgxy, Avgx2, Ver, Par, Diffy, width, Sig;
float Avg = 0, Diff = 0;

unsigned int pix(unsigned int x, unsigned int y) {
	return black[x * Width + y];
}

unsigned int nvis(unsigned int x, unsigned int y) {
	if((signed int)x < 0 || x >= Height || (signed int)y < 0 || y >= Width) return false;
	return light[x * Width + y];
}

unsigned char RGB_to_Black(rgb_t color) {
	return round((color.red + color.green + color.blue) / 3.0);
}

bool bfs(unsigned int x, unsigned int y) {
	std::queue <unsigned int> list;
	std::queue <unsigned int> list2;
	std::queue <unsigned int> list3;
	list.push(x * Width + y);
	light[x * Width + y] = false;
	unsigned int tp, t, i, j;
	Sum = 0, Sumx = 0, Sumy = 0, Sumxy = 0, Sumx2 = 0;
	float wgt;
	while(!list.empty()) {
		tp = list.front();
		list.pop();
		list2.push(tp);
		x = tp / Width;
		y = tp % Width;
		
		t = pix(x, y);
		wgt = t / 255.0;
		//printf("%d, %d\n", x, y);
		Sum += wgt;
		Sumx += wgt * x;
		Sumy += wgt * y;
		Sumxy += wgt * x * y;
		Sumx2 += wgt * x * x;
		
		if(nvis(x - 1, y)) list.push((x - 1) * Width + y), light[(x - 1) * Width + y] = false;
		if(nvis(x + 1, y)) list.push((x + 1) * Width + y), light[(x + 1) * Width + y] = false;
		if(nvis(x, y - 1)) list.push(x * Width + y - 1), light[x * Width + y - 1] = false;
		if(nvis(x, y + 1)) list.push(x * Width + y + 1), light[x * Width + y + 1] = false;
	}
	
	if(Avg * Size > Sum * LINEMUL) return false;
	
	//printf("Sum: %.3f, Sumx: %.3f, Sumy: %.3f, Sumxy: %.3f, Sumx2: %.3f\n", Sum, Sumx, Sumy, Sumxy, Sumx2);
	Avgx = Sumx / Sum, Avgy = Sumy / Sum, Avgxy = Sumxy / Sum, Avgx2 = Sumx2 / Sum;
	
	if(Sumxy - Avgx * Avgy * Sum == 0 || Sumx2 - Avgx * Avgx * Sum == 0) return false;
	
	Ver = (Sumxy - Avgx * Avgy * Sum) / (Sumx2 - Avgx * Avgx * Sum);
	Par = -1.0 / Ver;
	//printf("Avgx: %.3f, Avgy: %.3f, Avgxy: %.3f, Avgx2: %.3f, Ver: %.6f, Par: %.6f\n", Avgx, Avgy, Avgxy, Avgx2, Ver, Par);
	
	Diffy = 0;
	while(!list2.empty()) {
		tp = list2.front();
		list2.pop();
		list3.push(tp);
		x = tp / Width;
		y = tp % Width;
		wgt = pix(x, y) / 255.0;
		Diffy += wgt * abs(y - (Avgy + (x - Avgx) * Ver));
	}
	Diffy /= Sum;
	Diffy *= 2;
	width = Diffy / sqrt(Ver * Ver + 1);
	//printf("Diffy: %.6f Width: %.6f\n", Diffy, width);
	Sig = 0;
	while(!list3.empty()) {
		tp = list3.front();
		list3.pop();
		x = tp / Width;
		y = tp % Width;
		wgt = pix(x, y) / 255.0;
		Sig += abs(wgt * abs(y - (Avgy + (x - Avgx) * Ver)) - Diffy / 2);
	}
	Sig /= Sum;
	Sig = Sig / sqrt(Ver * Ver + 1);
	
	return true;
}

bool comp(point a, point b) {
	return a.Avgy < b.Avgy;
}

int main() {
	bitmap_image Input("RustOutput.bmp");
	freopen("linewidth.txt", "w", stdout);
	if(!Input) {
       fprintf(stderr, "Error - Failed to open: DivideInput.bmp\n");
       return 1;
    }
	// printf("Image read successfully!\n");
	
	// printf("Image size: %d * %d\n", Height = Input.height(), Width = Input.width());
	Height = Input.height();
	Width = Input.width();
	Size = Height * Width;
	
	unsigned int i, j, k, t;//i -> row, j -> column
	for(i = 0; i < Height; i++)
		for(j = 0; j < Width; j++) {
			black.push_back(t = RGB_to_Black(Input.get_pixel(j, i)));
			light.push_back(t > 0);
			cnt[t / CNTWIDTH]++;
			Avg += t;
		}
	Avg /= Size;
	for(i = 0; i < Height; i++)
		for(j = 0; j < Width; j++) {
			Diff += abs(pix(i, j) - Avg);
		}
	Diff /= Size;
	printf("灰度统计图：\n");
	for(i = 0; i < CNTNUM; i++) {
		printf("%03d - %03d: %.3f ", i * CNTWIDTH, (i + 1) * CNTWIDTH - 1, cnt[i] * 1.0 / Size);
		for(j = 0; j < cnt[i] * 100 / Size; j++) printf("*");
		printf("\n");
	}
	printf("Avg: %.6f, Diff: %.6f\n", Avg, Diff);
	
	int listnum = 0;
	point lines[std::max(Height, Width)];
	point pt;
	for(i = 0; i < Height; i++)
		for(j = 0; j < Width; j++) {
			if(nvis(i, j)) {
				if(!bfs(i, j)) continue;
				pt.Avgx = Avgx; pt.Avgy = Avgy; pt.Avgx2 = Avgx2; pt.Avgxy = Avgxy;
				pt.Sum = Sum; pt.Sumx = Sumx; pt.Sumy = Sumy; pt.Sumx2 = Sumx2; pt.Sumxy = Sumxy;
				pt.Par = Par; pt.Ver = Ver; pt.Diffy = Diffy; pt.width = width; pt.Sig = Sig;
				lines[listnum++] = pt;
			}
		}
	
	std::sort(lines, lines + listnum, comp);
	
	image_drawer draw(Input);
	draw.pen_width(1);
	for(int i = 0; i < listnum; i++) {
		pt = lines[i];
		Avgx = pt.Avgx; Avgy = pt.Avgy; Avgx2 = pt.Avgx2; Avgxy = pt.Avgxy;
		Sum = pt.Sum; Sumx = pt.Sumx; Sumy = pt.Sumy; Sumx2 = pt.Sumx2; Sumxy = pt.Sumxy;
		Par = pt.Par; Ver = pt.Ver; Diffy = pt.Diffy; width = pt.width; Sig = pt.Sig;
		
		printf("\nline %d:\n", i);
		printf("Sum: %.3f, Sumx: %.3f, Sumy: %.3f, Sumxy: %.3f, Sumx2: %.3f\n", Sum, Sumx, Sumy, Sumxy, Sumx2);
		printf("Avgx: %.3f, Avgy: %.3f, Avgxy: %.3f, Avgx2: %.3f, Ver: %.6f, Par: %.6f\n", Avgx, Avgy, Avgxy, Avgx2, Ver, Par);
		printf("Diffy: %.6f 粗细Width: %.6f 粗糙度Sig: %.6f\n", Diffy, width, Sig);
		printf("角度Angle: %.6f\n", atan(Par) / 3.14 * 180);
		
		draw.pen_color(255, 0, 0);
		draw.circle(Avgy, Avgx, 0);
		draw.circle(Avgy, Avgx, 1);
		draw.circle(Avgy, Avgx, 2);
		draw.circle(Avgy, Avgx, 3);
		
		draw.pen_color(0, 255, 0);
		draw.line_segment(Avgy, Avgx, Avgy + (Height - Avgx) * Ver, Height);
		draw.line_segment(Avgy, Avgx, Avgy - Avgx * Ver, 0);
		
		draw.pen_width(Sig);
		draw.pen_color(127, 0, 0);
		draw.line_segment(Avgy + Diffy + 1, Avgx, Avgy + (Height - Avgx) * Ver + Diffy + 1, Height);
		draw.line_segment(Avgy + Diffy + 1, Avgx, Avgy - Avgx * Ver + Diffy + 1, 0);
		draw.line_segment(Avgy - Diffy, Avgx, Avgy + (Height - Avgx) * Ver - Diffy, Height);
		draw.line_segment(Avgy - Diffy, Avgx, Avgy - Avgx * Ver - Diffy, 0);
		draw.pen_width(1);
		
		draw.pen_color(0, 0, 255);
		draw.line_segment(Avgy, Avgx, Avgy + width * Par / sqrt(1 + Par * Par), Avgx + width / sqrt(1 + Par * Par));
		draw.line_segment(Avgy, Avgx, Avgy - width * Par / sqrt(1 + Par * Par), Avgx - width / sqrt(1 + Par * Par));
	}
	
	//y - Avgy = Ver * (x - Avgx)
	//y = Ver * x - Ver * Avgx + Avgy
	//Ver * x - y + Avgy - Ver * Avgx = 0
	//x + Par * y - Par * Avgy - Avgx = 0
	//x = -Par * y + Par * Avgy + Avgx
	
	draw.pen_width(3);
	draw.pen_color(255, 255, 0);
	for(int i = 0; i < listnum - 1; i++) {
		point pt1 = lines[i], pt2 = lines[i + 1];
		Ver = (pt1.Ver + pt2.Ver) / 2;
		Par = -1.0 / Ver;
		float DisC = (pt2.Avgy - Ver * pt2.Avgx) - (pt1.Avgy - Ver * pt1.Avgx);
		float Dis = abs(DisC) / sqrt(Ver * Ver + 1);
		printf("\nline %d & %d:\n", i, i + 1);
		printf("Ver1: %.6f Ver2: %.6f\n", pt1.Ver, pt2.Ver);
		printf("Par: %.6f Ver: %.6f\n", Par, Ver);
		printf("DisC: %.6f 间距Dis: %.6f\n", DisC, Dis);
		printf("线宽Depression: %.6f\n", Dis + pt1.width + pt2.width);
		draw.line_segment(pt1.Avgy, pt1.Avgx, pt1.Avgy + Dis / sqrt(Ver * Ver + 1), pt1.Avgx - Ver * Dis / sqrt(Ver * Ver + 1));
		//Avgy1 - DisC = Par * x - Par * Avgx2 + Avgy2
		//Avgy1 - DisC + Par * Avgx2 - Avgy2 = Par * x
		//x = (Avgy1 - DisC - Avgy2) / Par + Avgx2
		
		//y = - Ver * Disc + Avgy
	}
	fclose(stdout);
	//printf("\n");
	// freopen("CON", "w", stdout);
	
	Input.save_image("DivideOutput.bmp");
	// printf("All things done!\nNow you can safely close the window!\n");
	//while(1);
	return 0;
}
